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Summary At the heart of any method for computational fluid dynamics lies 
the question of how the simulated fluid should be discretized. Traditionally, 
a fixed Eulerian mesh is often employed for this purpose, which in modern 
schemes may also be adaptively refined during a calculation. Particle-based 
methods on the other hand discretize the mass instead of the volume, yield- 
ing an approximately Lagrangian approach. It is also possible to achieve La- 
grangian behavior in mesh-based methods if the mesh is allowed to move with 
the flow. However, such approaches have often been fraught with substantial 
problems related to the development of irregularity in the mesh topology. Here 
we describe a novel scheme that eliminates these weaknesses. It is based on 
a moving unstructured mesh defined by the Voronoi tessellation of a set of 
discrete points. The mesh is used to solve the hyperbolic conservation laws of 
ideal hydrodynamics with a finite volume approach, based on a second-order 
Godunov scheme with an exact Riemann solver. A particularly powerful fea- 
ture of the approach is that the mesh-generating points can in principle be 
moved arbitrarily. If they are given the velocity of the local flow, a highly ac- 
curate Lagrangian formulation of continuum hydrodynamics is obtained that 
is free of mesh distortion problems, while it is at the same time fully Galilean- 
invariant, unlike ordinary Eulerian codes. We describe the formulation and 
implementation of our new Voronoi-based hydrodynamics, and we discuss a 
number of illustrative test problems that highlight its performance in practical 
applications. 



1 Introduction 

Numerical simulations have become an indispensable tool to study fluid dy- 
namics, especially in astrophysics where direct experiments are often impos- 
sible. However, it is not always clear whether the employed simulation algo- 
rithms are sufficiently accurate in real practical applications, and to which 
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extent numerical deficits may affect the final results. It therefore remains an 
important task to continue to critically test the numerical methods that are 
in use, and to develop new approaches with the goal to reach better accuracy 
at comparable or even lower computational cost. 

In astrophysics, a variety of fundamentally quite different numerical meth- 
ods for hydrodynamical simulations are i n use, the most prominent ones are 
smoothed particle hydrodynamics (SPH ; iLucvl . Il977 : Gingold fc Monagharl 
19771 : iMonagharll 19921: Springel l2010bh and Eulerian mesh -based hydrody- 
namics (e.g. iTorol . Il997l : iLeVeauel . 120021 ; IStone et all l2008h with (optional) 
adaptive mesh refinement (AMR) . A particular challenge in astronomy is the 
need to calculate self-gravitating flows, which often tend to cluster strongly 
under gravity, producing a huge dynamic range in density and length scales 
that can only be treated efficiently with spatially adaptive resolution. An 
important reason for the popularity of SPH lies in the fact that such an adap- 
tivity is automatically built into this method, whereas achieving it in adaptive 
mesh refinement codes requires substantial effort. 

It has become clear over recent years that both SPH and AMR suffer 
from fundamental problems that make them inaccurate in certain regimes. 
Indeed, these methods sometimes yield conflicting results even for basic cal- 
culations that only consider non-radiative hydrodynamics (e.g. iFrenk et al 



19991 ; lAgertz et all 120071 ; iTasker et all 120081 ; iMitchell et ail . |2009j). SPH codes 



have comparatively poor shock resolution, offer only low-order accuracy for 
the tr eatment of contact discontinuities, and suffer from subsonic velocity 
noise (|Abel I 201 ll ) Worse, they appea r to suppress fluid instabilities under 
certain conditions (jAgertz et al.l . 120071 ) , as a result of a spurious surface ten- 
sion and inaccurate gradient estimates across density jumps. On the other 
hand, Eulerian codes are not free of fundamental problems either. They do 
not produce Galilean-invariant results, w hich can make their accuracy sensi- 
tive t o the presence of bulk velocities (e.g. Wadslev et al. . 2008 : Tasker et al.l . 
2008). Another concern lies in the mixing inherent in multi-dimensional Eu- 



lerian hydrodynamics. This provides for an implicit sourc e of entropy, with 



sometimes unclear consequences (e.g. Wadslev et al. . 2008f ). 



There is hence substantial motivation to search for new hydrodynamical 
methods that improve on these weaknesses of the SPH and AMR techniques. 
In particular, we would like to retain the accuracy of mesh-based hydrody- 
namical methods (for which decades of experience have been accumulated in 
computational fluid dynamics) , while at the same time we would like to outfit 
them with the Galilean-invariance and geometric flexibility that is character- 
istic of SPH. The principal idea for achieving such a synthesis is to allow the 
mesh to move with t he flow itself. This is an obvious and old idea fBraun fc 
Sambridge"Tl995l; IGnedinl. Il995l: IWhitehurstl. Il995 l; iMavriplisl . Il997l; IXul. Il997 : 
Hassan et all Il998l ; iPenl . 1 19981 ; iTrac fc Perl |200J) , but one fraught with many 
practical difficulties that have so far prevented widespread use of any of the 
few past attempts to int roduce moving -mes h me t hods in astrophysics and 
cosmology. For example, IGnedinl ( 199a ) and IPenl ( 19981 ) presented moving- 
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Fig. 1. Example of a Voronoi and Delaunay tessellation in 2D, with periodic bound- 
ary conditions. The red circles show the generating points of the Voronoi tessellation, 
which is drawn with solid lines. Its topological dual, the Delaunay triangulation, is 
overlaid with thin dashed lines. 

mesh hydrodynamic algorithms which relied on the continuous deformation 
of a Cartesian grid. However, the need to limit the maximum allowed grid dis- 
tortions severely impacts the flexibility of these codes for situations in which 
the mesh becomes heavily distorted, and special measures are required to 
let the codes evolve cosmological density fields into a highly clustered state. 
In general, mesh tangling (manifested in 'bow-tie' cells and hourglass like 
mesh motions) is the traditional problem of such attempts to simulate multi- 
dimensional hydrodynamics in a Lagrangian fashion. 

In this contribution, we describe a new formulation of continuum hydro- 
dynamics based on an unstructured mesh. The mesh is defined as the Voronoi 
tessellation of a set of discrete mesh-generating points, which are in principle 
allowed to move freely. For the given set of points, the Voronoi tessellation 
of space consists of non-overlapping cells around each of the sites such that 
each cell contains the region of space closer to it than to any of the other 
sites. Closely related to the Voronoi tessellation is the Delaunay tessellation, 
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the topological dual of the Voronoi diagram. Both constructions have already 
been widely used for natu ral neighbor interpolatio n and geometric analysis 
of cosmic structures (e.g. Ivan de Wevgaertl. 119941: ISambridge et all . Il995l ; 
Schaap &: van de Wevgaertl. l2000l: IPelupessv et al. . 20031; van de Wevgaert & 



Schaap, 20091 ). In 2D. the Delaunav tessellation for a given set of points is a 
triangulation of the plane, where the points serve as vertices of the triangles. 
The defining property of the Delaunay triangulation is that each circumcircle 
around one of the triangles of the tessellation is not allowed to contain any 
of the other mesh-generating points in its interior. This empty circumcircle 
property distinguishes the Delaunay triangulation from the many other trian- 
gulations of the plane that are possible for the point set, and in fact uniquely 
determines the triangulation for points in general position. Similarly, in three 
dimensions, the Delaunay tessellation is formed by tetrahedra that are not 
allowed to contain any of the points inside their circumspheres. 

As an example, Figure [T] shows the Delaunay and Voronoi tessellations for 
a small set of points in 2D, enclosed in a box with imposed periodic boundary 
conditions. The midpoints of the circumcircles around each Delaunay trian- 
gle form the vertices of the Voronoi cells, and for each line in the Delaunay 
diagram, there is an orthogonal face in the Voronoi tessellation. 

The Voronoi cells can be used as control volumes for a finite-volume for- 
mulation of hydrodynamics, using the same principal ideas for reconstruction, 
evolution and averaging (REA) steps that are commonly employed in many 
Eulerian techniques. However, as we will see it is possible to consistently in- 
clude the mesh motion in the formulation of the numerical steps, allowing the 
REA-scheme to become Galilean-invariant. Even more importantly, due to the 
mathematical properties of the Voronoi tessellation, the mesh continuously de- 
forms and changes its topology as a result of the point motion, without ever 
leading to the dreaded mesh-tangling effects that are the curse of traditional 
moving mesh methods. We note that the approach we describe here is quite 
different from attempts to formulate fluid particl e models based on Vor onoi 
cells (e.g. iHietel et all l2000l : ISerrano et all 120051 ; IrlefJ fc Springel . |2010| ). or 
mesh- free finite volume approaches ( Junkl . 12002 ). The former are similar in 
spirit to SPH and typically maintain a constant mass per particle, whereas 
our scheme is really closely related to ordinary mesh codes - except that the 
mesh is fully dynamic. 

With illustrative test problems we shall later show that the resulting for- 
mulation of hydrodynamics performs rather well on a number of test problems, 
featuring very high accuracy in the treatment of shocks, shear waves, and fluid 
instabilities. In particular, it can give better results than fixed-mesh Eulerian 
hydrodynamics, thanks to much reduced advection errors. It also offers much 
higher accuracy than SPH when an equal number of particles/cells is used, 
making it highly attractive as a possible alternative to currently employed 
SPH and AMR schemes in astrophysics 

This article is structured as follows. In Section [21 we formulate continuum 
hydrodynamics on the Voronoi mesh, based on a finite-volume ansatz and 
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a second-order accurate extension of Godunov's method. In Section [3J we 
briefly discuss time integration and implementation aspects. We then turn to 
a discussion of a number of basic hydrodynamical tests in Section SI chosen to 
highlight some of the principal advantages and properties of the new approach. 
We note that a more extensive discussion of code tests and of the algorithmic 
i mplementation of the new scheme in the parallel AREPO code can be found 
in ISpringell Finally, we summarize and discuss our main findings in 

Section [SJ 



2 A finite volume discretization of the Euler equations 
on a moving Voronoi mesh 

The Euler equations are conservation laws for mass, momentum and energy 
that take the form of a system of hyperbolic partial differential equation. They 
can be written in compact form by introducing a state vector 




U = | pv | = | pv j (1) 

pu + \pv 2 

for the fluid, where p is the mass density, v is the velocity field, and e = 
u + v 2 /2 is the total energy per unit mass, u gives the thermal energy per 
unit mass, which for an ideal gas is fully determined by the temperature. 
These fluid quantities are functions of the spatial coordinates x and time t, 
i.e. U = U(x, t), but for simplicity we will typically refrain from explicitly 
stating this dependence in our notation. Based on U, we can define a flux 
function 




F(U)= | u / | (2i 
with an equation of state 

P = ( 7 - l)pu (3) 

that gives the pressure of the fluid. The Euler equations can then be written 
in the compact form 

<9U 

-+V-F.0, (4) 

which emphasizes their character as conservation laws for mass, momentum 
and energy. 

Over the past decades, a large variety of different numerical approaches 
to s olve this coupled set of part ial differential equations have been developed 
(see Torol . 1997 ; LeVeauel . 2002 , for comprehensive expositions). We will here 



employ a so-called finite-volume strategy, in which the discretization is carried 
out in terms of a subdivision of the system's volume into a finite number of 
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disjoint cells. The fluid's state is described by the cell-averages of the conserved 
quantities for these cells. In particular, integrating the fluid over the volume 
Vi of cell i, we can define the total mass m^, momentum pi and energy 
contained in the cell as follows, 



Qi= Pi = / UdV. (5) 




With the help of the Euler equations, we can calculate the rate of change of 
Qi in time. Converting the volume integral over the flux divergence into a 
surface integral over the cell results in 



dQi 
df 



/ [F(U) - Uw T ] dn. (6) 

JdVi 



Here n is an outward normal vector of the cell surface, and w is the velocity 
with which each point of the boundary of the cell moves. In Eulerian codes, 
the mesh is taken to be static, so that w = 0, while in a fully Lagrangian 
approach, the surface would move at every point with the local flow velocity, 
i.e. w = v. In this case, the right hand side of equation (j6|) formally simplifies, 
because then the first component of Qi, the mass, stays fixed for each cell. Un- 
fortunately, it is normally not possible to follow the distortions of the shapes 
of fluid volumes exactly in multi-dimensional flows for a reasonably long time, 
or in other words, one cannot guarantee the condition w = v over the entire 
surface. In this case, one needs to use the general formula of equation |6j), as 
we will do in this work. 

The cells of our finite volume discretization are polyhedra with flat polyg- 
onal faces (or lines in 2D). Let Ay describe the oriented area of the face 
between cells i and j (pointing from i to j). Then we can define the averaged 
flux across the face i-j as 

■•Y. I [F(U) - Uw T ] dAy, (7) 

and the Euler equations in finite-volume form become 

3 

We obtain a manifestly conservative time discretization of this equation by 
writing it as 

Q(n+l) = Q (»)_ zit ^ A .F(™+ 1 / 2 ) ) (g) 
3 

where the Fy are now an appropriately time-averaged approximation to the 
true flux Fy across the cell face. The notation Q| is meant to describe the 
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Fig. 2. Sketch of a Voronoi mesh and the relevant geometric quantities that enter 
the flux calculation across a face. In (a), we show the the mesh-generating points i\ 
and Tj of two cells i and j. The face between these two cells has a center-of-mass 
vector fij, which in general will be offset from the mid-point rriij of the two points. 
In (b), we illustrate the two velocity vectors w, and Wj associated with the mesh- 
generating points. These are normally chosen equal to the gas velocity in the cells, 
but other choices are allowed too. The motion of the mesh-generating points uniquely 
determines the motion of the face between the cells. Only the normal velocity w is 
however needed for the flux computation in the rotated frame x' , y . 



state of the system at step n. Note that Fy = — Fjj, i.e. the discretization is 
manifestly conservative. 

Evidently, a crucial step lies in obtaining a numerical estimate of the fluxes 
Fij , and a good fraction of the literature on computational fluid dynamics is 
concerned with this problem. This issue is particularly important since the 
most straightforward (and perhaps naive) approach for estimating the fluxes, 
namely simply approximating them as the average of the left and right cell- 
centered fluxes catastrophically fails and invariably leads to severe numerical 
integration instabilities that render such a scheme completely useless in prac- 
tice. 

Many modern schemes for estimating the fluxes in a stable fashion are 
descendants of Godunov's method, which revolutionized the field. By solv- 
ing an exact or approximate Riemann problem at cell boundaries, Godunov's 
method allows the correct identification of the eigenstructure of the local so- 
lution and of the upwind direction, which is crucial for numerical stability. 
While Godunov's original method offers only first order accuracy and is rela- 
tively diffusive, it can be extended to higher-order accuracy relatively simply, 
and in many different ways. We employ the MUSCL-Hancoc k scheme (V an 
Leer. Il984j ; iTorol . Il997t Ivan Leerl . I 2006T ) . which is a well-known and relatively 
simple approach for obtaining second-order accuracy in space and time. This 
scheme is also popular in astronomy and used in several state-of-the art Eule- 



rian c odes fe.g. lFromang et al.l . l2006l ; lMignone et all 120071 ; ICunningham et al 



2009) . In its basic form, the MUSCL-Hancock scheme involves a slope- limited 



piece- wise linear reconstruction step within each cell, a first order prediction 
step for the evolution over half a timestep, and finally a Riemann solver to 
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estimate the time-averaged inter-cell fluxes for the timestep. After the fluxes 
have been applied to each cell, a new averaged state of the cells is constructed. 
This sequence of steps in a timestep hence follows the general REA approach. 

Figure [5] gives a sketch of the geometry involved in estimating the flux 
across the face between two Voronoi cells. Truly multidimensional Riem ann 
solv ers have been developed recently ( Wendrofit 19991: Brio et all 200ll: Bal- 



, 2010), but it is unclear whether they can be readily adapted to our 



sara, 

complicated face geometry. We therefore follow the common approach and 
calculate the flux for each face separately, treating it as an effectively one- 
dimensional problem. Since we do not work with Cartesian meshes, we cannot 
use operating splitting to deal with the individual spatial dimensions. Rather 
we use a method where all the fluxes are computed in one step, and are then 
collectively applied to calculate the change of the conserved quanti ties in a 



cell. This unsplit approach implicitly accounts for "corner fluxes" ([Colellal . 
1990) needed to recover second-order accuracy through the half-step predic- 
tion step in our MUSCL-Hancock scheme (which exploits the primitive form 
of the Euler equations). For defining the Riemann problem normal to a cell 
face, we rotate the fluid state into a suitable coordinate system with the x'- 
axis normal to the cell face (see sketch) . This defines the left and right states 
across the face, w hich we pass to an exact Riemann solver. The latter is imple- 



mented following iTorol (Il997h with an extension to treat vacuum states, but 



it could easily be substituted with an approximate Riemann solver for higher 
performance, if desired. We note that in multi dimensions the transverse ve- 
locities are also required in the Riemann problem in order to identify the 
correct upwind transverse velocity, which is important for an accurate treat- 
ment of shear. Once the flux has been calculated with the Riemann solver, we 
transform it back to the lab frame. 

A further important point concerns the treatment of the allowed motion 
of cell surfaces in our scheme. In order to obtain stable upwind behavior, 
the Riemann problem needs to be solved in the frame of the moving face. 
This is important as the one-dimensional Riemann problem is not Galilean- 
invariant in the following sense: Suppose left and right state at an interface 
are described by (pl, Pl,vi,) and (pr,Pr,vr), for which the Riemann solver 
returns an interface state (pp, Pf,vf) that is the basis for the flux estimate. For 
example, the mass flux across the interface is then given by pfVf- Consider now 
a velocity boost v applied both to the left and the right side. The new Riemann 
problem is given by (pl,-Pl, + v) and (pr,-Pr,wr + v), and will return a 
flux estimate p' F v' F . However, in general this will yield p' F v F ^ Pf(vf + v), 
which implies that the calculated flux vector is not Galilean invariant. 

In our new hydrodynamical scheme, each timestep involves the following 
basic steps: 

1. Calculate a new Voronoi tessellation based on the current coordinates 
of the mesh generating points. This also gives the centers-of-mass Sj of 
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each cell, their volumes Vi , as well as the areas Ay and centers fy of all 
faces between cells. 

2. Based on the vector of conserved fluid variables associated with each 
cell, calculate the 'primitive' fluid variables Wj = (pi, vy Pi) for each cell. 

3. Estimate the gradients of the density, of each of the velocity components, 
and of the pressure in each cell, and apply a slope-limiting procedure to 
avoid overshoots and the introduction of new extrema. 

4. Assign velocities w t to the mesh generating points. 

5. Evaluate the Courant criterion and determine a suitable timestep size At. 

6. For each Voronoi face, compute the flux Fy across it by first determining 
the left and right states at the midpoint of the face by linear extrapolation 
from the cell midpoints, and by predicting these states forward in time 
by half a timestep. Solve the Riemann problem in a rotated frame that is 
moving with the speed of the face, and transform the result back into the 
lab- frame. 

7. For each cell, update its conserved quantities with the total flux over its 
surface multiplied by the timestep, using equation ([9]). This yields the new 
state vectors Q^™ +1 ^ of the conserved variables at the end of the timestep. 

8. Move the mesh-generating points with their assigned velocities for this 



For the sake of definiteness, we now briefly describe the most important details 
of these different steps. 

2.1 Gradient estimation and linear reconstruction 

According to the Green-Gauss theorem, the surface integral of a scalar func- 
tion over a closed volume is equal to its gradient integrated over the same 
volume, i.e. 



This suggests one possible way to estimate the mean gradient in a Voronoi 
cell, in the form 



where 0(fy) is the value of <fi at the centroid fy of the face shared by cells i 
and j, and Ay is a vector normal to the face (from j to i), with length equal 
to the face's area. Based on the further approximation 



timestep. 




(10) 




(11) 



(12) 



this provides an estimate for the local gradient. Note that with the use of 
equation (fl"2"j) . the gradient of cell i only depends on the values <j)j of neigh- 
boring cells, but not on itself. While the estimate (jlll) can be quite generally 
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applied to arbitrary tessellations, due to the use of only one Gauss point per 
face it is also relatively inaccurate and is not exact to linear order in general. 

For the special case of Voronoi cells, it is however possible to obtain a 
considerably better gradient estimate with little additional effort. The key is 
to carr y out the surface integral more accurately. It can be shown ( Serrano & 
Espaiiol" T200ll : ISpringell . l2010al ) that the gradient estimate 



1 , / r , , 1 Cii <Pi+ ' 



W>« - - ( - 4h] ^ - ^ ^ ) (13) 

is exact to linear order, independent of the locations of the mesh-generating 
points of the Voronoi tessellati on. Here w e followed the notation of Serrano 
& Espanol ( 200ll ) and defined c,j as the vector from the midpoint between 



i and j to the center-of-mass of the face between i and j. Without the term 
involving Cjj this gradient estimate is the same as the simpler Green-Gauss 
estimate. However, retaining this extra term leads to significantly better accu- 
racy, because the gradient estimate becomes exact to linear order for arbitrary 
Voronoi meshes. In practice, we shall therefore always use this gradient es- 
timation in our MUSCL-Hancock scheme for the Euler equations, where we 
calculate in this way gradients for the 5 primitive variables (p,v Xl v y ,v z , P) 
that characterize each cell. 

The result (TT31 has also a n interesting relation to the formulae obtained 
by ISerrano fc: Espanol ( 2001 ) for the partial derivatives of the volume of a 



Voronoi cell with respect to the location of one of the po ints. As Serrano & 



Espanol (|200lf ) have shown, the derivative of the volume of a Voronoi cell due 



to the motion of a surrounding point is given by 



Furthermore, they show that 



3 



(15) 



Using these relations, and noting that according to the Gauss theorem we 
have 

#£^= ' (16) 

because the summation is just the surface integral of a constant function, we 
can also write the estimate for the gradient of <fi at more compactly as 

3 

An interesting corollary of the above is that provided <j>(r) varies only linearly, 
the sum 
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S = Y i 4>(Ti)V i (18) 

i 

approximates the integral J 4>(r) dV exactly, independent of the positions of 
the points that generate the Voronoi tessellation. 

In our approach, we use the gradients estimated with equation (fT3|) for a 
linear reconstruction in each cell around the center-of-mass. For example, the 
density at any point r G Vi of a cell is estimated as 

p(r) = Pi + {Vp) i ■ (r - Si ), (19) 

where Sj is the center of mass of the cell. Note that independent of the mag- 
nitude of the gradient and the geometry of the Voronoi cell, this linear re- 
construction is conservative, i.e. the total mass in the cell rrii is identical to 
the volume integral over the reconstruction, m, = J v p(r)d 3 r. An alternative 
choice for the reference point is to choose the mesh-generating point instead 
of This is the more natural choice if the cell values are known to sample the 
values of the underlying field at the location of the mesh-generating points, 
then the reconstruction is exact to linear order. However, our input quanti- 
ties are cell-averages, which correspond to linear order to the values of the 
underlying field sampled at the center-of-masses of the cells. For this reason 
we prefer the center-of-mass of a cell as reference point for the reconstruction. 

Nevertheless, this highlights that large spatial offsets between the center- 
of-mass of a cell and its mesh-generating point are a source of errors in the 
linear reconstruction. It is therefore desirably to use "regular" meshes if pos- 
sible, where the mesh-generating points lie close to the center-of-mass; such 
meshes minimize the errors in the gradient estimation and the linear recon- 
struction. Or in other words, we would like our Voronoi meshes to be relatively 
close to so-called centroidal Voronoi meshes, where the mesh-generating points 
lie exactly in the center of mass of each cell. As we discuss in bit more detail 
later, we have developed an efficient method for steering the mesh motion 
such that this regularity condition can be approximately maintained at all 
times. 

2.2 Slope limiting procedure 

In smooth parts of the flow, the above reconstruction is second-order accurate. 
However, in order to avoid numerical instabilities the order of the reconstruc- 
tion must be reduced near fluid discontinuities, such that the introduction of 
new extrema by over- or undershoots in the extrapolation is avoided. This 
is generally achieved by applying slope limiters that reduce the size of the 
gradients near local extrema, or by flux limiters that replace the high-order 
flux with a lower order version if there are steep gradients in the upstream 
region of the flow. 

We here generalize the original MUSCL approach to an unstructured grid 
by enforcing monotonicity with a slope limiting of the gradients. To this end 
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we require that the linearly reconstructed quantities on face centroids do not 



excee d the maxima or minima among all neighboring cells (jBarth fc Jesperson 



1989). Mathematically, we replace the gradient with a slope-limited gradient 

(V0>; = on (V<^ , (20) 

where the slope limiter < on < 1 for each cell is computed as 

OLi = mm(l,<if>ij). (21) 

Here the minimum is taken with respect to all cells j that are neighbors of 
cell i, and the quantity ipij is defined as 

f (C ax - <k)IMu for Mij > 

i>a = { (C in - <k)IMa for Mij < (22) 

[ 1 for Afcj = 

where A<j>%j — (V0)^ • (fy — Sj) is the estimated change between the centroid 
fij and the center of cell i, and 0™ ax = max(0j) and ^™ ln = max((4j) are the 
maximum and minimum values occurring for cj> among all neighboring cells 
of cell i, including % itself. We note that this slope limiting scheme does not 
strictly enforce the total variation diminishing property, which means that 
(usually reasonably small) post-shock oscillations can sometimes still occur. 



2.3 Setting the velocities of the mesh generators 

A particular strength of the scheme we propose here is that it can be used 
both as an Eulerian code, and as a Lagrangian scheme. The difference lies only 
in the motion of the mesh-generation points. If the mesh-generating points are 
arranged on a Cartesian mesh and zero velocities are adopted for them, our 
method is identical to a second-order accurate Eulerian code on a structured 
grid. Of course, one can equally well choose a different layout of the points, in 
which case we effectively obtain an Eulerian code on an unstructured mesh. 
The real advantage of the new code can be realized when we allow the mesh 
to move, with a velocity that is tied to the local fluid speed. In this case, we 
obtain a Lagrangian hydrodynamical code, which has some unique and impor- 
tant advantages relative to an Eulerian treatment. It is however also possible 
to prescribe the mesh motion through an external flow field, for example in or- 
der to smoothly concentrate resolution towards particular regions of a mesh, or 
to realize rotating meshes. Unlike other arbitrary Lagrangian-Eulerian (ALE) 
fluid dynamical methods, the method proposed here does however not rely 
on remapping techniques to recover from distortions of the mesh once they 
become severe, simply because the Voronoi tessellation produced by the con- 
tinuous motion of the mesh-generating points yields a mesh geometry and 
topology that itself changes continuously in time, without mesh-tangling ef- 
fects. 
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The most simple and basic approach for specifying the motion of the mesh 
generators is to use 



i.e. the points are moved with the fluid speed of their cell. This Lagrangian 
ansatz is clearly appropriate for pure advection and in smooth parts of the 
flow. However, in this scheme there is no mechanism built in that tries to 
improve the regularity of the Voronoi mesh in case the mean mass per cell 
should develop substantial scatter around a desired mean value, or if cells 
with high aspect ratios occur. If desired, such tendencies of a growing mesh 
irregularity can be counteracted by adding corrective velocity components 
to the primary mesh velocities of equation (f2"3")l . There are many different 
possibilities for how exactly to do this, and we consider this freedom a strength 
of the formalism. In Section[3T2j we will discuss a simple regularization method 
that we have found to be very effective. 

2.4 Flux computation 

An important aspect of our approach is that the specified velocities of the 
mesh-generating points fully determine the motion of the whole Voronoi mesh, 
including, in particular, the velocities of the centroids of cell faces (see sketch 
in Fig. [2J . This allows us to calculate the Riemann problem in the rest- frame 
of each of the faces. 

Consider one of the faces in the tessellation and call the fluid states in the 
two adjacent cells the 'left' and 'right' states. We first need to determine the 
velocity w of the face based on the velocities wl and wr of the two mesh- 
generating points associated with the face (they are connected by a Delaunay 
edge). It is clear that w has a primary contribution from the mean velocity 
(wl + wr)/2 of the points, but there is also a secondary contribution w' from 
the residual motion of the two points relative to their center of mass. This 
residual motion is given by w R = — = (wr — wl)/2, and we need to 
determine its impact on the motion of the face centroid. The components of 
w R and Wl parallel to the line connecting the centroid f of the face with the 
midpoint m of the two mesh-generating points Tl and Tr induce a rotation 
of the face around the point m. We are only interested in the normal velocity 
component of this motion at the centroid of the face. This can be easily 
computed as 



(23) 



, _ (w L - w R ) ■ [f - (tr + r L )/2] (r R - r L ) 
|rR-r L | |r R -r L |' 

The full velocity w of the face is then given by 



(24) 



Wr -I- W L 



+ w'. 



(25) 



w 



2 



We now calculate the flux across the face using the MUSCL-Hancock ap- 
proach, with the important difference that we shall carry out the calculation 
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in the rest-frame of the face. It is convenient to do this in the primitive vari- 
ables (p,v,P), where we first transform the lab-frame velocities of the two 
cells to the rest-frame of the face by subtracting w, 

WU = W L , R - w . (26) 



We then linearly predict the states on both side to the centroid of the face, 
and also predict them forward in time by half a timestep. This produces the 
states 



w" w. 



L,R - VV L,R T- —q^- 



(f-s L , R )+ — 

L,R Ub 



f- (27) 

L,R Z 



The spatial derivatives dW /dr are known, and given by the (slope-limited) 
gradients of the primitive variables that are estimated as described in Sec- 
tion 12.11 Note that the gradients are unaffected by the change of rest-frame 
described by Eqn. (|26[) . The partial time derivate dW/dt can be replaced by 
spatial derivatives as well, based on the Euler equations in primitive variables, 
which are given by 

dw aw 

_ +A( W)_ = 0, (28) 

where A is the matrix 



A(W) = v l/p . (29) 




Having finally obtained the states left and right of the interface, we need to 
turn them into a coordinate system aligned with the face, such that we can 
solve an effectively one-dimensional Riemann problem. The required rotation 
matrix A for the states only affects the velocity components, viz. 



= A W" R = A 3D W" , (30) 




where A3D is an ordinary rotation of the coordinate system, such that the 
new x-axis is parallel to the normal vector of the face, pointing from the left 
to the right state. 

With these final states, we now solve the Riemann problem, and sample 
the self-similar solution along x/t = 0. This can be written as 

W F = i? ic mann(W^,W R '), (31) 

where -Riemann is a one-dimensional Riemann solver, which returns a solution 
for the state of the fluid Wp on the face in primitive variables. We now 
transform this back to the lab-frame, reversing the steps above, 
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Wi ab = vi ab = A-'W F + w . (32) 



Finally, we can use this state to calculate the fluxes in the conserved variables 
across the face. Here we need to take into account that the face is moving 
with velocity w, meaning that the appropriate flux vector in the lab frame is 
given by 

/ KVlab - W) \ 

F = F(U)-Uw T = /9v lab (v lab -w) T + P , (33) 

\peiab(viab - w) + Pv lab / 

where U is the state Wi a b expressed in the conserved variables, and ei a b = 
v i 2 a b/2 + flab/ [(7 — l)piab]- The scalar product of this flux vector with the 
normal vector of the face gives the net flux of mass, momentum, and energy 
that the two adjacent, moving cells exchange. It is the flux of equation (|33l) 
that can finally be used in the conservative updates of each cell, as described 
by equation ([9]). 

For the above formulation, it is straightforward to show that the changes 
in the conserved quantities in a cell are Galilean invariant; any Galilean boost 
is effectively simply absorbed into the motion of the face. 



3 Time integration and implementation aspects 
3.1 Time integration 

For hydrodynamics with a global timestcp, we employ a simplified CFL 
timestep criterion in the form 

= Ccfl R \ ,, (34) 

to determine the maximum allowed timestep for a cell i. Here Ri is the ef- 
fective radius of the cell, calculated as Ri = (3V^/47r) 1 / 3 from the volume 
of a cell (or as Ri = (Vi/n) 1 / 2 from the area in 2D), under the simplifying 
assumption that the cell is spherical. The latter is normally a good approxi- 
mation, because we steer the mesh motion such that the cell-generating point 
lies close to the center-of-mass of the cell, which gives it a "roundish" polyhe- 
dral shape. Ccfl < 1 is the Courant-Friedrichs-Levy coefficient (usually we 
choose Ccfl — 0.4 — 0.8), Cj = \/jP/p is the sound speed in the cell, and 
l v 'il = l v » — w i\ 1S the velocity of the gas relative to the motion of the grid. 
In the Lagrangian mode of the scheme, the velocity |v^| is close to zero and 
usually negligible against the sound speed, which means that larger timesteps 
than in an Eulerian treatment are possible, especially if there are large bulk 
velocities in the system. 
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If the code is operated with a global timestep, we determine the next 
system timestep as the minimum 

At = min AU (35) 

i 

of the timestep limits of all particles. 

It is also possible to implement an individual timestep scheme, where the 
timestep conditions of different cells are treated in a more flexible fashion. 
This can greatly improve the computational efficiency in many applications. 
For example, in cosmological simulations, a large dynamic range in densities 
quickly occurs as a result of gravitational clustering. Accordingly, the local 
dynamical times can vary by orders of magnitude. It has therefore long be- 
come common practice to use individual timesteps for the collisionless N-body 
problem, a te chnique that has also been extended to hy drodynamical SPH sim- 



problem, a technique that has also been extended to hydrodynamical brti sim- 
ulations (e.g. Katz et al. . 19961 : Springel et al. L l200lh . We have implemented 



such a method also for the moving-mesh scheme, based on a discretization 
of the allowed timestep sizes into a power-of-two hierarchy. Unlike the ap- 
proach taken in AMR simulations, where rehned grid patches are typically as 
a whole subcycled in time by a constant factor, we impose no such restriction 
on our mesh, i.e. in principle each cell can be evolved with its own timestep, 
constrained only to the power-of-two hierarchy of allowed timestep sizes. To 
maintain a fully conservative character of the scheme, we evolve each face with 
the smaller timestep of the two adjacent c ells. Fu ll details of this individual 
timestep scheme can be found in Springel! ( 2010a ). 



3.2 Mesh regularity 

As seen in Figure [1] Voronoi meshes may sometimes look quite "irregular" , 
in the sense that there is a significant spread in sizes and aspect ratios of the 
cells, especially for disordered point distributions. While this is not a problem 
of principle for our approach, it is clear that the computational efficiency will 
normally be optimized if regions of similar gas properties are represented with 
cells of comparable size. Having a mixture of cells of greatly different volumes 
to represent a gas of constant density will restrict the size of the timestep 
unnecessarily (which is determined by the smallest cells), without giving any 
benefit in spatial resolution (which will be limited by the largest cells in the 
region). 

As we have seen, it is also desirable to have cells where the center-of- 
mass lies close to the mesh-generating point, because this minimizes errors in 
the linear reconstruction and limits the rate at which mesh faces turn their 
orientation during mesh motion. Below, we will discuss one possible approach 
for steering the mesh motion during the dynamical evolution such that, if 
desired, mesh regularity in the above sense can be achieved and m aintained. 

In so-called centroidal Voronoi tessellations ( Okabe et al. . 200fj| ). the mesh- 



generating points coincide with the center-of-mass of all cells. There is an 
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amazingly simple algorithm known as Lloyd's method ( Llovdl . [l982) to obtain 



a centroidal Voronoi tessellation starting from an arbitrary tessellation. One 
simply moves the mesh-generating points of the current Voronoi tessellation 
to the center-of-masses of their cells, and then reconstructs the Voronoi tessel- 
lation. The process is repeated iteratively, and with each iteration, the mesh 
relaxes more towards a configuration in which the Voronoi cells appear quite 
'round' (in 2D they form a honeycomb-like mesh) and have similar volume - 
a centroidal Voronoi tessellation. 

Inspired by this algorithm, we employ a simple scheme to improve, if 
needed, the local shape of the Voronoi tessellation during the dynamical evo- 
lution. We simply augment equation (I23[) with an additional velocity com- 
ponent, which is designed to move a given mesh-generating point towards 
the center-of-mass of its cell. There are different possibilities to parameterize 
such a corrective velocity. One approach that we found to work quite well in 
practice is to add a correction velocity whenever the mesh-generating point 
is further away from the center-of-mass of a cell than a given threshold, ir- 
respective of the actual velocity field of the gas. To this end, we associate a 
radius Ri = (SVi/Air) 1 / 3 with a cell based on its volume (or area in 2D). If 
the distance di between the cell's center-of-mass Sj and its mesh-generating 
point Yi exceeds some fraction 77 of the cell radius Ri, we add a corrective 
term proportional to the local sound speed Cj of the cell to the velocity of 
the mesh-generating point. This effectively applies one Lloyd iteration (or a 
fraction of it) to the cell by repositioning the mesh-generating point onto the 
current center-of-mass, ignoring other components of the mesh motion. In or- 
der to soften the transition between no correction and the full correction, we 
parameterize the velocity as 



but the detailed width of this transition is unimportant. In very cold flows 
the sound speed may be so low that the correction becomes ineffective. As an 
alternative, we therefore also implemented an option in our code that allows 
a replacement of c s (si — Yi)/di in equation (|36|) with (s$ — Yi)/ At. This more 
aggressive approach to ensure round cells generally works very well too, but 
has the disadvantage to depend on the timestepping. Our typical choice for 
the threshold of the correction is r] — 0.25, and we usually set % = 1A 
i.e. the correction is, if present, applied in full over the course of one timestep. 
Smaller values of 77 can be used to enforce round cell shapes more aggressively, 
if desired. 

The above scheme is usually quite effective in maintaining low aspect ratios 
and a regular mesh at all times during the evolution. However, the criterion 
for detecting cells that should get a correction velocity is not triggered if a 
mesh is strongly stretched or compressed in one direction; then the centers of 




(36) 
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mass of cells can still be close to their mesh-g enerating points, but the aspect 
ratio of cells can be very high. We have found ([Vogelsberger et all 1201 ll ) that a 
simple alternative criterion treats such situations much better. To this end we 
determine for each cell the maximum angle under which any of the faces of the 
cell is seen from its mesh-generating point. If this angle lies above a prescribed 
threshold value, the mesh-correction component to the velocity is invoked, just 
as above. This approach will effectively try to prevent that a mesh-generating 
point gets too close to an outer wall of a cell, which simultaneously ensures 
that the displacement from the center-of-mass and the aspect ratio stay small. 
We have also found that with this criterion the mesh-correction motions are 
required more rarely, hence we have made this our default choice for general 
simulations with the moving mesh approach. In any case, it is important to 
note that the correction velocities are still Galilean-invariant, and they vanish 
most of the time, so that the mesh-generating points will usually be strictly 
advected with the local fluid velocity. 

We point out that there is an important difference of thi s approach com- 
pared with the mesh regularization technique presented in Hefi fc Springell 
(2010). \n the finite volume approach discussed here, one may in principle 
move the mesh-generating points in nearly arbitrary ways. Maintaining a 
good mesh is therefore comparatively str aightforward, as descri bed above. 
In contrast, the Voronoi particle model of lHefi fc Springel ( 2010h dictates a 
particular equation-of-motion for the mesh-generating points, where one is 
not allowe d to simply add so me mes h correction ve locities. As a way out. Hefi 
& Springel d201^" suggested to modify the underlying Lagrangian in a tricky 
way in order to automatically build in corrective motions into the dynamics 
of the mesh-generating points, but this approach is not equally flexible as the 
one we can use here. 



3.3 Implementation aspects 

The scheme described thus fa r has been impl emented in the AREPO code, 
which is described in detail in ISpringell (|2010ah . A central aspect of the code 



is a fast engine for the generation of Delaunay and Voronoi meshes. To this 
end an incremental insertion algorithm is used both in 2D and 3D, which 
also allows partial mesh constructions, as needed in our individual timestep 
approach. For reasons of memory and run-time efficiency, we have written 
our own low-level mesh-construction routines instead of using a library such 
as CGAL. The mesh construction is parallelized for distributed memory ma- 
chines. We use a spatial domain decomposition in which each domain is 
mapped to a single processor, which then first constructs its part of the mesh 
independently of the other CPUs, and then exchanges and inserts additional 
'ghost' particles as needed to make sure that the Voronoi cells of all local par- 
ticles are complete, i.e. that their geometry is identical to the one that would 
be found in a fiducial global mesh constructed in serial. 
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In order to robustly treat degenerate cases (for example when more than 
three points lie on a common circle), we work with a computational vol- 
ume that is mapped to double precision numbers in the interval [1,2[. For 
IEEE arithmetic, the mantissa of these numbers effectively defines a one-to- 
one mapping of all representable floating point numbers in this range to the 
space of 53-bit integers. We then evaluate geometric predicates with fast ordi- 
nary double precision arithmetic, but always monitor the maximum round-off 
error. If the outcome of a predicate is not guaranteed to be correct as a result 
of round-off errors, we compute the predicate exactly with long-integer arith- 
metic based on the mantissas corresponding to the floating point numbers. 
We find that this approach is both robust and still quite fast. 

Finally, we would like to mention that our moving-mesh approach can 
quite easily be coupled to self-gravity using similar algorithms as are often 
employed in particle-based SPH codes. In fact, in the AREPO code we use a 
simil ar TreePM solver for gravity as employed in the GADGET code ([Springe! 
20051 ). In our approach, the Voronoi cells are treated effectively as point masses 



with a gravitational softening length set equal to the fiducial radius of the cell, 
as estimated from its volume. The tree-code has the advantage of being highly 
efficient also for strongly clustered particle configurations, and it can be easily 
adapted to individual timcstcp integration. 



4 Illustrative test problems 

We now discuss a number of simple test problems that show the performance 
of the moving mesh approach and illustrate its specific strengths. In some 
cases, we will compare directly to SPH simulations based on the same initial 
conditions. Also, we discuss differences in the solutions when the mesh is 
instead kept stationary, in which case our method behaves equivalently to a 
standard Eulerian scheme with second-order accuracy in space and time. 

4.1 Riemann and Sod-shock problems 

Arguably the most important basic test problems of hydrodynamical codes 
consist of one-dimensional Riemann problems. In the Riemann problem, two 
piece-wise constant states, each characterized by density, pressure and ve- 
locity, are brought into contact with each other, and their subsequent time 
evolution is then followed. If the initial velocities are zero, one deals with the 
special case of a Sod-shock problem. The Riemann problem can be solved 
analytically for an ideal gas, and the solution consists of a set of three self- 
similar waves that emerge from the initial discontinuity. There is in general 
one contact wave in the middle, sandwiched on either side by either a shock 
or a rarefaction. The ability of a hydrodynamical method to accurately treat 
different Riemann problems is fundamental for the ability of the scheme to 
capture complex hydrodynamical phenomena. 
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Fig. 3. Different one-dimensional Riemann problems, calculated with a resolution 
of 100 points in the unit domain, for a gas with adiabatic index 7 = 1.4. The three 
columns show results for the initial conditions of the problems 1, 2 and 3 as specified 
in the text. Symbols represent the hydrodynamical quantities of the Voronoi cells, 
while the solid lines give the analytic solutions for density, velocity and pressure, 
from top to bottom. 



In order to highlight the ability of our moving Voronoi mesh to accu- 
rately represent Riemann problems we consider in Fig. [3] the results for three 
Riemann problems, simulated in ID with 100 initially equally spaced points 
in the unit domain. For problem one, the initial conditions are given by 
(px,Pi,Vi,p 2 ,P2,v 2 ) = (1.0,1.0,1.0,0.125,0.1,0.0), for problem 2 the cor- 
responding values are (1.0, 0.4, —2.0, 1.0, 0.4, 2.0), and for problem 3 they are 
(1.0, 1000.0, 0.0, 1.0, 0. 01, 0.0). The se values describe the same problems as dis- 
cussed in the book bv lTorol(il997l ). and correspond to a moderate Sod shock 
tube, a strong double rarefaction, and a very strong shock. 
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As we can see in Fig. [3] from the comparison to the analytic solution, 
the moving mesh approach captures the solutions of these Riemann problems 
rather accurately. The contact discontinuities and shocks are quite sharp, 
with only negligible post-shock oscillations. The only significant error occurs 
in the nearly evacuated region of the strong rarefaction, where the simulated 
temperature is too high. However, this is a common error of Eulerian codes 
when applied to this problem. We also see that the spatial resolution varies 
at the end, since the points have moved with the flow. In particular, the 
resolution has become quite low in the low density region that develops in the 
middle of the double rarefaction, while it has increased on the right hand side 
of the contact discontinuities in the two Sod-shock problems. We note that the 
accuracy with which the analytic solution is recov ered is considera bly better 
than with SPH for the same initial conditions (see Springel . l2010bh . 



4.2 Isentropic vortex convergence 

We next turn to a test of a non-trivial multi-di mensional fluid problem with a 



smoo th solution, the isentropic vortex problem ()Yee et al.l . l2000l ; ICalder et al 



2002). This problem is particularly useful for verifying whether our method 
does indeed show second-order convergence, despite the presence of strong de- 
formations and topological changes of the mesh and the use of small velocity 
components to keep the mesh nice and regular, as described above. Previ- 
ously, second-order convergence of the AREPO cod e has only been explicitly 
demonstrated for ID sound waves (ISpringeilillOah . which is a comparatively 
simple problem where no mesh twisting occurs. The 2D Gresho vortex test on 
the other hand ( Gresho fc Chan . 19901 : Liska fc Wendrofi . 2003 ) shows only a 
convergence rate of —1.4 as a function of the number of cells per dimension 
( Springe i l2010al ). both in our moving mesh code and other fixed mesh codes 



that have second-order accuracy, like ATHENA (jStone et all 120081 ) . However 



this can be understood as a result of the presence of discontinuities in the 
Gresho problem. 

Yee et al. describe the setup of a perfectly smooth vortex, which has 

an analytic, time-invariant solution. To realize this problem, we adopt a box 
of extension [—5, 5] 2 in 2D, with periodic boundaries everywhere. The initial 
distribution of mesh-generating points is adopted as a regular Cartesian grid. 
The velocity field is specified as 



/? f l-r 2 
v x {x,y) = -V^ CX P I — 2 — 

/3 f l-r 2 
v y {x,y) = x—cic-p I - 



(37) 
(38) 



where r 2 = x 2 + y 2 . The density and thermal energy per unit mass are calcu- 
lated from 

T(x, y) ee P/p = Too — {l ~ y exp (l - r 2 ) (39) 
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Fig. 4. L2 error norm at time t = 8.0 as a function of resolution for the isentropic 
vortex test (left panel). The dashed line is a power law with L2 oc iV -2 . The panels 
on the right show the mesh in the region [—4, 4] x [—4, 4] around the center at times 
t = 1.5 (top) and t — 8.0 (bottom), for the run at resolution 80 x 80. 



as p — T 1 ^ 7-1 ) and u = T/(<y — 1). For these choices, the entropy P/p 1 is 
exactly constant everywhere, and the solution is time- independent. We adopt 
7 = 1.4, a vortex strength f3 — 5.0, and = 1. We realize the initial con- 
ditions by integ rating over th e fields in each grid cell to obtain the conserved 
variables, as in ICalder et al.l ([20021 ) . We then calculate the evolution of the 
vortex with our moving mesh code for different resolutions until time t — 8.0, 
at which point the vortex has rotated more than once, and the mesh has been 
thoroughly sheared in the region of the vortex. 

In Figure IH we consider the L2-norm of the numerically obtained density 
field at the final time relative to the analytic solution, as a function of reso- 
lution. We use N 2 = 40 2 , 80 2 , 160 2 , 320 2 , 640 2 , and 1280 2 initial mesh cells. 
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Reassuringly, the error declines accurately as a powerlaw, with L2 oc iV~ 2 , 
which is the expected convergence rate for a se cond-order ac c urate scheme. 
This convergence rate has also been reached by ICalder et all (|2002h for the 
FLASH code, but unlike for this code, the error in our approach is completely 
independent on whether or not the vortex has an additional bulk velocity. 
We note that our result also disagrees with the conjecture that an additional 
(second) mesh-construction per time step would be needed to reach th is con- 
vergence rate for multi-dimensional flow (jDuffell &: MacFadvenl . l201lh . 



4.3 Rayleigh- Taylor instabilities 

In multi-dimensional flows, a further class of important hydrodynamical phe- 
nomena besides acoustic waves and the non-linear waves related to Riemann 
problems (shocks, contact discontinuities and rarefaction waves) appears. 
These are so-called fluid instabilities, such as the Rayleigh- Taylor or Kelvin- 
Helmholtz instabilities. They are highly important for producing turbulence 
and for inducing mixing processes between different phases of fluids. 

The Rayleigh- Taylor instability can arise in stratified layers of gas in an 
external gravitational field. If higher density gas lies on top of low-density gas, 
the stratification is unstable to buoyancy forces, and characteristic finger-like 
perturbations grow that will mix the fluids with time. To illustrate this insta- 
bility and simultaneously show the motion of the mesh in our Voronoi based 
code, we illustrate in Figure [5] the evolution of a single Rayleigh- Taylor mode, 
calculated at the deliberately low resolution of 12 x 36 cells. The simulation 
domain is two-dimensional, with extension [0.5, 1.5] and periodic boundaries 
at the vertical boundaries, and solid walls at the bottom and top. There is an 
external gravitational field with acceleration g = —0.1, and the bottom and 
top halves of the box are filled with gas of density p = 1 and p = 2, respec- 
tively. The gravitational forces are balanced by an initial hydrostatic pressure 
profile of the form P(y) = Pa + {y — 0.75) g p(y) with Pq and 7 = 1.4. To seed 
the perturbation, one mode is excited with a small velocity perturbation of 
the form v y (x,y) — wq[1 — cos(47ra)][l — cos(47ry/3)], where wo — 0.0025. 

As can be clearly seen in the time evolution shown in Figure [5j the 
Rayleigh- Taylor instability is captured well by the moving-mesh method even 
at this low resolution. What is particularly interesting is that the sharp bound- 
ary between the phases can be maintained for relatively long time during the 
early evolution of the instability, simply because the contact discontinuity is 
not smeared out as it bends, thanks to the mesh's ability to follow this motion 
in an approximately Lagrangian fashion. A Eulerian approach with a station- 
ary mesh on the other hand would automatically wash out the boundary due 
to advection errors, involving some spurious mixing of the fluids. 

This fundamental improvement of the moving mesh code with respect to 
a fixed mesh approach becomes clearer in Figure |H1 Here we compare a high- 
resolution version of the Rayleigh- Taylor instability between the moving-mesh 
approach and the same calculation carried out with a stationary Cartesian 
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Fig. 5. Rayleigh- Taylor instability calculated at low resolution with the moving- 
mesh approach. A denser fluid lies above a less dense fluid in an external gravitational 
field. The hydrostatic equilibrium of the initial state is buoyantly instable. The three 
frames show the time evolution of the density field of the system at times t = 5.0, 
10.0, and 15.0, after a single mode has been perturbed to trigger the stability, as 
described in the text. 

mesh. Here 1024 x 1024 cells have been used in the unit domain, [—0.5, 0.5] 2 , 
and the instability was triggered by adding small random noise to the y- 
velocity field, of the form v y (x, y) = A [1 + cos(27ry)]/2, where A is a random 
number in the interval [—0.05,0.05]. While the instability shows a similar 
overall growth rate in both cases, eventually leading to full turbulence in the 
box, there are also striking differences. Whereas the calculation with the fixed 
mesh produces a lot of intermediate density values due to the strong mixing 
of the phases on small scales, the moving mesh approach maintains finely 
stratified regions where different layers of the fluid phases have been folded 
over each other. The contact discontinuities between these layers can be kept 
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t = 2.5, moving mesh 


t = 7.5, moving mesh 


t = 12.5, moving mesh 


1 = 2.5, fixed mesh 


t= 7.5, fixed mesh 


1 = 12.5, fixed mesh 



Fig. 6. Rayleigh- Taylor instability calculated at high resolution with 1024 x 1024 
points in the unit domain. The instability is here seeded by small random noise 
added to the velocity field. The top and bottom rows compare the time evolution 
for calculations with a moving and a stationary mesh, respectively. 



sharp by the code even when they are moving relative to the rest-frame of 
the box. We think that this behavior is much more faithful to the underlying 
hydrodynamical flow. In the early phase of the growth, it also appears as if 
small-scale RT fingers grow somewhat too quickly in the Eulerian case as a 
result of grid alignment effects. 

4.4 Kelvin-Helmholtz instabilities 



Another important fluid instability arises in shear flows, the so-called Kelvin- 
Helmholtz (KH) instability. Whenever there is a discontinuity in the shear 
velocity across a fluid interface, wave-like transverse perturbations across the 
interface will grow into characteristic wave-like patterns. This instability is 
ultimately behind the generation of waves on lakes and oceans when wind 
streams over the water. The KH instability is ubiquitous in complex flows 
and plays a prominent role in the generation of turbulence. 

It has recently been found that the simulation of KH instabilities can 
be quite problematic in SPH, with the growth b eing suppressed wh en the 



density jump across the fluid discontinuity is large (|Agertz et all 120071 ). This 



has triggered a flurry of activity in the recent literature on SPH, trying to 
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Fig. 7. Growth rate of the KH instability for two different density ratios, 1 : 2 (left 
panel) or 1 : 10 (right panel) . The red lines show the results for the moving mesh code 
when initial conditions with a sharp density jump are used, either with an initial 
perturbation amplitude of vo = 0.001 or vq — 0.01. The dashed blue lines are the 
corresponding results for a stationary mesh. The dotted lines give the exponential 
growth expected from linear perturbation theory. Finally, the thin red and thin blue 
dashed lines are the results obtained when the initial discontinuity is washed out in 
the initial conditions. Time is given in units of the KH growth timescale tkh- 



improve on this behavi or dPricel. 120081 iRead et all l2010l : iHefi fc Springe! 
20101 : lJunk et all l201fj| : lAbej . l201lh . 

We here show a basic KH test in a two-dimensional set-up, comparing 
our moving-mesh approach against the traditional fixed-mesh approach. For 
dcfinitcness, we fill a box with periodic boundaries and unit length on a side 
with gas of density P2 — 2 in the horizontal middle stripe, and the rest with 
density pi = 1. The middle region is moving to the right with velocity v x = 0.5, 
the rest of the gas moves to the left with velocity v x = —0.5. The initial 
distribution of the mesh-generation points is a Cartesian grid of resolution 
256 x 256. We seed an initial perturbation by adding an additional component 



v y {x,y) = vq sin(fc:r) 



(40) 



to the velocity field, where vo is a small number. We choose k = 2 x (2w/L). 
Hence the Fourier spectrum of the v y field contains in the beginning only the 
k x = 2 mode. The growth of this mode for t > can then be conveniently 
measured through Fourier transforms of the velocity field. 

In the left panel of Figure a number of measurements of the numerical 
KH growth rate for a density ratio of 1 : 2 are summarized. The red lines show 
the result for our new moving-mesh method; the upper line is for an initial 
perturbation amplitude of vq — 0.01, while the lower line is for vq = 0.001. 
The blue dashed lines correspond in both cases to a fixed Cartesian mesh of 
the same resolution. The linear theory growth rate v y oc exp(i/rKH) is shown 
as dotted lines, where 
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moving mesh, sharp contact, t = 1.0 



fixed mesh, sharp contact, t = 1.0 





I fixed mesh, softened contact, t = 2.0 




Fig. 8. Kelvin Helmholtz instability computed with different initial conditions, and 
for a moving or a fixed mesh. In the panels on top, the initial contact discontinuity 
was sharp between adjacent cells. In contrast, in the bottom row it was smoothed 
out. 



tkh = | 1 r . (41 

\V2 - v 1 \k v /pip 2 

is the KH growth timescale for an inviscid gas. Because the density is initially 
not perturbed self-consistently with the velocity field, it takes first a bit of time 
before the instability develops, but then the moving-mesh solution follows the 
expected linear theory growth rate quite nicely for a while. Eventually the 
growth slows down as the mode saturates and the non-linear evolution of the 
KH instability ensues. 

In the fixed-mesh case, the results are somewhat less clean and depend 
more strongly on the initial perturbation amplitude. Visual inspection of den- 
sity maps during the time evolution reveals that not only the excited mode 
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starts to grow but also shorter wavelength modes. This is not too surprising 
since small wavelength perturbations grow fastest in the Kelvin- Helmholtz in- 
stability, and representing the initial sharp density jump implicitly involves a 
spectrum of small waves. The latter are prone to outgrow the larger-scale seed 
perturbation once the discontinuity starts to be misaligned with the principal 
coordinate axes, in which case additional small-scale perturbations are seeded 
at mesh corners. 

This effect can be repaired if the initial c ontact discontinuity is washed 
out, as advocated by iRobertson et ah ( 2010| ). Then also in the fixed- mesh 
case only the excited mode grows and a more stable result is obtained. The 
latter is shown as thin dashed blue line in Figure [3 However, in this case one 
does not reach the full growth rate that is expected analytically for the (sharp) 
instability with this wavenumber. If one also applies the same smoothing to 
the initial conditions of the moving-mesh run, one obtains essentially the same 
result for the growth rate of the exited mode, which is shown with a thin red 
line. As soon as the initial discontinuity is smooth enough to be resolved by 
several mesh cells, it hence appears as if it would not make a difference whether 
one uses a moving or a fixed-mesh. This is however not true. If one waits long 
enough it is seen that the moving-mesh code resolves secondary KH billows for 
which the fixed-mesh approach appears to be already too diffusive. This can 
be seen in the density maps of Figure|8l where the bottom two panels compare 
the density field of the KH test at time t = 2.0 for the moving and the fixed 
mesh approaches, using smoothed initial conditions. The top two panels on 
the other hand give the same comparison when a sharp initial discontinuity 
is used instead. In the latter case, it is clearly seen that more 'wrong' modes 
grow in the simulation with a stationary mesh, because here a misalignment 
of the the sharp boundary with the mesh triggers larger seed perturbations 
on small scales than for the moving mesh. The question whe ther this initial 
condition is somehow 'allowed' or not ( Robertson et al. . 2010l ) is moot in our 



view. Both codes are started from identical initial conditions, and hence the 
comparison tests how susceptible the codes are to the growth of numerically 
seeded small-scale perturbations in this situation. 

Finally, we have also considered the same test with a density jump of I : 10, 
which was realized by raising the density of the middle stripe to p2 — 10. The 
pressure was increased by a factor of 5 to P — 12.5 in order to ensure that 
the flow stays supersonic and the Mach number of the shear flow is not much 
changed. The right panel in Figure [7] shows the corresponding results for the 
growth rate. Qualitatively, the results closely follow those obtained for the 
1 : 2 density ratio, except that the fixed-mesh results show a markedly too 
fast overall growth rate in this case. 



4.5 Shock-cloud interaction 



We finally consider two problems that involve the interaction of strong shocks 
with fluid instabilities, which is important in many astrophysical applications. 
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Fig. 9. Comparison of the time evolution of the cloud disruption test in 2D with 
the moving mesh code (left column) and with SPH (right column). The panels show 
the density field at diff erent times, as labeled. The same set-up and the identical 
initial conditions as in ISpringell (|2005l ) have been used. In the t = 4.0 frame of the 
moving-mesh calculation, a small rectangle marks a region that is shown enlarged 
at the bottom left, with the Voronoi-mesh overlaid. 



First, we repeat a test presented in the GADGET code paper ([Springe! 120051 ) 
as an advanced test of SPH. Here a strong shock wave of Mach number 10 
strikes an initially overdense cloud with density p = 5 that is embedded at 
pressure equilibrium in a tenuous hot phase of density p — 1 and pressure 
P = 1 (with 7 = 5/3). Figure |9] shows the time evolution of the system in 2D, 
comparing the moving mesh results (left row) with the SPH result obtained 
with the GADGET code for the identical initial conditions (right row). As 
the shock strikes the cloud, it is compressed and accelerated. A complicated 
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system of reflected and interacting shocks develops, and in the flow around 
the cloud, vortices are generated by the baroclinic term. With time, these 
vortices tend to at least partially disrupt the cloud. 

In comparing the moving-mesh and the SPH results a number of interesting 
observations can be made. First, the density field in the smooth regions is 
noticeably noisy in the SPH calculation when compared with the moving- 
mesh approach. Also, the shock waves are not as sharp and crisp as in the 
Voronoi-based code, even though the global flow features are clearly very 
similar in both cases. Arguably the most important difference is however that 
the cloud is shredded much more in the moving-mesh simulation, while the 
SPH result shows a large degree of coherence of the cloud debris. In fact, little 
"droplets" of dense gas remain that are eventually advected downstream in 
the SPH calculation, showing no tendency to mix further with the background 
gas. This is presumably related to a spurious surface tension effect in SPH 
across contact discontinuities with large density jumps. 

In the bottom left panel of the time-sequence shown in Figure [9j we have 
marked a small region with a black rectangle. In order to illustrate the geom- 
etry of the Voronoi mesh in this simulation, this region is shown enlarged at 
the bottom of Fig. HI with the mesh overlaid. It can be seen that the higher 
density region in the top right is populated with smaller Voronoi cells than the 
lower density region at the bottom left, as a result of the Lagrangian character 
of the scheme. 

Finally, we turn to a related simulation problem in 3D. whic h has become 
known as the 'blob-test'. First carried out in lAgertz et al. ( 20071 ). this consists 



of a three-dimensional overdense sphere that is put into a low-density back- 
ground gas that streams su personically with re spect to the cloud. We adopt 



the same parameters as in lAgertz et al.l ([20071 ). and use the original mesh- 



based initial conditions of this test as made available on the interne10, at three 
different resolutions equal to 32x32x64, 64x64x128, and 128x128x256. Sim- 
ilar to the two-dimensional shock-cloud interaction problem discussed above, 
the supersonic head wind leads to the development of a shear-flow over the 
surfa ce of the cloud, which produces disrupting KH instabilities. In Agertz 
et al. (l2007f l it was found that the SPH calculations would only lead to an 
incomplete destruction of the cloud, while the considered Eulerian mesh-code 
predicted a complete destruction of the cloud after a relatively short time. 
The latter was measured in terms of the mass-fraction of the original cloud 
that was still denser than 0.64 times the density of the cloud in the initial 
conditions and colder than 0.9 times the temperature of the background gas. 

In Figure 1101 we show our results for the remaining cloud mass fraction 
as a function of time for three different resolutions, both for the moving- 
mesh code and for the equivalent calculation with a stationary mesh. We find 
that the moving-mesh approach appears to give converged results already 
at lower resolution than the Eulerian approach. It also consistently shows a 



1 They can be downloaded at http://www.astrosim.net 
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Fig. 10. The remaining mass fraction of a dense gas blob as a function of time 
when it is pu t into a supersonic head wind, corresponding to the three-dimensional 
'blob test' of lAgertz et alJ l|2007t ). The red solid, dashed, and dotted lines show our 
results for the moving Voronoi-mesh, with initial resolutions of 128 x 128 x 256, 
64 x 64 x 128, and 32 x 32 x 64 cells, respectively. The blue lines give the results 
when the mesh is kept stationary instead. 



slightly higher residual cloud mass fraction than the fixed mesh calculation. 
Both effects can be understood as a result of the Galilean-invariance and the 
considerably lower advection errors of the moving mesh code. However, it is 
clear that both methods are qualitatively consistent and predict a complete 
disruption of the cloud after a timescale of t ~ 3tkh- Because SPH gives a 
lower mass loss at la t e times and does not produce a complete disruption of the 
cloud (|Agertz et all 120071 ; iHefi fc Springe! l2010h , this reinforces the concern 
that gas stripping out of dense system can be systematically underestimated 
in SPH. 



5 Discussion 

We have described a novel hydrodynamical scheme on an unstructured mesh 
that is constructed as the Voronoi tessellation of a finite set of mesh-generating 
points. The points are free to move during the time evolution, allowing the 
mesh to seamlessly follow the flow and to change its spatial resolution fully 
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adaptively. Thanks to the mathematical properties of the Voronoi tessellation 
there are no mesh-tangling or mesh-twisting effects since the motion of the 
mesh-generating points induces a continuous deformation of the mesh, without 
producing topological artefacts. 

Our approach represents a finite volume discretization of the Euler equa- 
tions on an unstructured mesh with second-order accuracy both in space and 
time, without the need to invoke an artificial viscosity. Unlike ordinary Eule- 
rian codes, the new method is fully Galilean-invariant, which is a very substan- 
tial advantage especially in simulations with large bulk flows. In particular, 
this property implies high accuracy for contact discontinuities and leads to a 
substantial reduction of advection errors when compared to traditional Eule- 
rian schemes. Indeed, these advantages of the moving-mesh approach can be 
readily demonstrated with test problems involving fluid instabilities or mov- 
ing shock waves. The Voronoi-based approach also avoids preferred spatial 
directions and offers flexibility in the treatment of boundary conditions. For 
example, curved boundaries or moving interfaces can be readily implemented. 

In this paper, we have not included a discussion of the technical aspects 
of implementing the scheme in t he new p a rallel c osmological code AREPO, 
which is described in full detail in ISpringel ( 2010a ). This implementation has 
already been applied to ti mely problems in cosmolog ical structure formation, 



such as ga laxy format i on (IVqgelsberger et all 120111 ) or the formation of the 



first stars ( Gr eif et al 



as radiative transfer 



namics ( Pakmor et al 



2011a 



13) . Recently , impo rtant physics extensions such 



Pet kwafc Springel l201ll ) and ideal magnetohydrody- 
201 if ) have been implemented in AREPO as well. The 



code is hence becoming an increasingly powerful alternative to more estab- 
lished simulation techniques in astrophysics. 

This is despite the considerable computational cost that the Voronoi mesh 
construction entails, and despite the complicated bookkeeping code that is 
required for the mesh management in parallel. In 3D hydrodynamics, our 
Voronoi code at present is about a factor of 2 slower for the same number of 
resolution elements than a SPH code (if 64 smoothing neighbours are used). 
Compared to a Eulerian fixed mesh code, the speed difference is about a fac- 
tor 3-4 (part of this difference also stems from the about twice larger average 
number of faces for our polyhedral cells compared with cubical cells in the 
Cartesian case) . A further discussion of the speed differen ce and the scalabil- 



ity of t he code for large cosmological applications is given in lVogelsberger et al 
( 20111 ). We note however that once self-gravity is added, the relative speed 
differences are much reduced, as often a sufficiently accurate calculation of 
gravity over a large dynamic range is more expensive than the hydrodynam- 
ics itself. Further note that our new technique reaches a given accuracy in 
many problems already at a lower resolution than SPH and fixed mesh codes, 
outweighing its higher complexity and making it hence also attractive from 
the point of view of computational efficiency. 
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